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Abstract-We describe a simple and highly efficient and accurate radiative transfer technique for 
computing bidirectional reflectance of a macroscopically flat scattering layer composed ot 
nonabsorbing or weakly absorbing, arbitrarily shaped, randomly oriented and randomly 
distributed particles. The layer is assumed to be homogeneous and optically semi-infinite, and 
the bidirectional reflection function (BRF) is found by a simple iterative solution of the 
Ambartsumian’s nonlinear integral equation. As an exact solution of the radiative transfer 
equation, the reflection function thus obtained fully obeys the fundamental physical laws of 
energy conservation and reciprocity. Since this technique bypasses the computation of the 
internal radiation field, it is by far the fastest numerical approach available and can be used as an 
ideal input for Monte Carlo procedures calculating BRFs of scattering layers with 
macroscopically rough surfaces. Although the effects of packing density and coherent 
backscattering are currently neglected, they can also be incorporated. The FORTRAN 
implementation of the technique is available on the World Wide Web at 
http://www.giss.nasa.gov/~crmim/brf.html and can be applied to a wide range of remote sensing, 
engineering, and biophysical problems. We also examine the potential effect of ice crystal shape 
on the bidirectional reflectance of fiat snow surfaces and the applicability of the Henyey- 
Greenstein phase function and the 5-Eddington approximation in calculations for soil surfaces. 
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1. INTRODUCTION 

Many remote sensing, engineering, and biophysical applications rely on accurate knowledge of 
the bidirectional reflection function (BRF) of layers composed of discrete, randomly positioned 
scattering particles (e.g., Refs. 1-46). Theoretical computations of BRFs for plane-parallel 
particulate layers are usually reduced to solving the radiative transfer equation (RTE) using one 
of existing exact or approximate techniques. Since some semi-empirical approximate 
approaches such as the Hapke model 13 are notorious for their low accuracy, crude violation of 
the energy conservation law, and ability to produce unphysical results, 28 ’ 31 the use of numerically 
exact solutions of RTE has gained justified popularity. For example, the computation of BRFs 
for particulate layers with macroscopically flat surfaces in Refs. 5, 17, and 19-22 is based on the 
adding-doubling technique. 47 ' 48 while Refs. 9 and 10 employ the discrete ordinate method. 49 
BRF computations for layers with undulated (macroscopically rough) surfaces are more 
complicated and often may have to rely on time-consuming Monte Carlo procedures. This 
approach is especially inefficient for optically thick, weakly absorbing media (e.g., snow and 
desert surfaces at visible wavelengths) since a photon may undergo many internal scattering 
events before it exits the medium or is absorbed. However, particulate layers with undulated 
surfaces can often be represented as collections of locally flat tilted facets characterized by the 
BRF found from the traditional plane-parallel RTE. In this way the Monte Carlo procedure 
could be used only to evaluate the effects of surface shadowing and multiple surface reflections, 
thereby bypassing the time-consuming ray tracing inside the medium and providing a great 
saving of CPU time. 

A further saving of computer resources can be achieved by using a more efficient 
technique for solving the plane-parallel RTE for a semi-infinite medium than the 
adding/doubling and discrete ordinate methods. Since many natural and artificial particulate 
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layers can be considered optically semi-infinite and homogeneous, one can find the BRf directly 
by solving the Ambartsumian's nonlinear integral equation' 11 using a simple iterative 
technique.' 1 ' 3 * In this way, the computation of the internal radiation field is avoided (cf. Rets. 
47—49) and the computer code becomes highly efficient and very accurate and compact. In the 
following sections, we discuss in detail numerical aspects and the computer implementation of 
this technique, examine the applicability of the Henyey-Greenstein phase function and the 5- 
Eddington approximation in BRF and flux calculations for soil surfaces, and describe sample 
applications demonstrating the potential effect of ice crystal shape on the bidirectional 
reflectance of flat snow surfaces. The last section summarizes the results of the paper and 
outlines further potential improvements of the model. 

2. COMPUTATIONAL TECHNIQUE 

We assume that the scattering layer is optically semi-infinite, has a macroscopically flat 
surface, and is composed of randomly distributed and randomly oriented particles of arbitrary 
shape. For simplicity, we ignore polarization effects and use intensity as the only physical 
characteristic of light. This so-called scalar approximation is not necessarily good for Rayleigh 
scattering, 53 ’ 54 but appears to be sufficiently accurate for particles with sizes comparable to and 
larger than the wavelength. 53 To describe the geometry of light scattering, we use a right-handed 
spherical coordinate system with the z axis directed along the outward normal to the surface (Fig. 
1). The direction of light propagation is specified by the couple (»,< p), where u = -cosS, 3 e 
[0°,180°] is the zenith angle, and (p € [0°,360°] is the azimuth angle. The azimuth angle is 
measured in the clock-wise direction when looking in the positive z direction. Note that u < 0 for 
upwelling radiation and u > 0 for downwelling radiation. We also define (.1 = \u\. The surface is 
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illuminated by a beam of unpolarized light incident in the direction (ito, <po - 0). The intensity of 
the retleeted radiation is given by 

/{-(.L,(p) = |.i 0 /?(n4i () .cp)f, (I) 

where /?(j.i,po. ( p) is the bidirectional reflection function and nF is the incident flux per unit area 
perpendicular to the incident beam. We ignore the effects of packing density, coherent 
backscattering, and shadow hiding (see Section 5) and find the reflection function as a 
numerically exact solution of the conventional radiative transfer equation. 50 ' 36 Specifically, we 
expand in a Fourier series in azimuth, 

ft(p,u 0 ,cp) = /? 0 (ji,n 0 ) + 2^7? m (p,p 0 )cosw(p , (2) 

OT = t 

and solve numerically the Ambartsumian’s nonlinear integral equation 50,5 ' 

(p + p 0 )/? m (p,p 0 ) = ^i 5 OT (-p,p 0 ) 

l 

ai> 0 )dp' 

2 0 

1 (3) 

+ lpji?'”(p,p')p m (p',p 0 )dp' 

0 

l 1 

0 0 

where gj is the single scattering albedo and / y "(|.i,p') are Fourier components of the phase 
function: 

F(p,p',(p-(p') = / >0 (|i,p') + 2^/ > '”(p,p , )cosm(cp-(p')- (4) 

m-l 
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The upper summation limit in Eq. (2) is chosen such that the absolute numerical accuracy of the 


Fourier expansion of the reflection function is better than a small predefined number e (e.g., e — 
I0 -4 ). 

The Fourier components of the phase function are given by 

p" W) = (-ir z N a,/>;; 0 (u)0(fO, (5) 

s-m 

where /^„(. r ) are generalized spherical functions 5 ' ' s closely related to associate Legendre 
functions (Section 3.2) and a s are expansion coefficients appearing in the standard expansion of 
the phase function P(Q ) in Legendre polynomials P s (x) = Pq 0 (x) : 


•^max , 

/>(©) = I a,/>(cos®), c< 0 = 1 , (6) 

s - 0 

where 0 is the scattering angle and s max is chosen such that all expansion coefficients with 
s > s max are smaller than O.le. Note that wz max <s max and we assume the following standard 
normalization of the phase function: 


- j>(0)sin0d© = l . (7) 

2 o 

If the expansion coefficients a 5 are known, then one can easily compute the Fourier components 

of the phase function via Eq. (5) and finally solve Eq. (3) using the method of simple iterations. 

It has been found that the method of simple iterations works very well for all m > 0. 
Furthermore, convergence is reasonably fast for m — 0 provided that the particles are absorbing 
(55 < 1). However, iterations converge very slowly or even may diverge for nonabsorbing or 
weakly absorbing particles (1 - 55 < 1). 51 It has been proved mathematically that this behavior is 
explained by the non-uniqueness of solutions of Eq. (3). :i9 To ameliorate this convergence 


5 



problem, DIugach and Yanovitskij ?l suggested to modify the /?"'(p,p () ) value after each iteration 
bv e n fo rein” the so-called Sobolev-van de Hulst relation 


I 

/(-Li) = 2 1 /f°(n,|x 0 )/(|.i 0 )|J. 0 dp 0 . 


( 8 ) 


The function /(//) is the solution of the equation 


f i 

/(»)( 1 - = — JV )^°(“>»') , 


(9) 


-i 


in which the so-called diffusion exponent k is found by satisfying the normalization condition 


rn 

1 


T I 

y J/(») du = l. 


( 10 ) 


After Eq. (3) is solved for each m, the evaluation of the Fourier series of Eq. (2) finalizes 
the process of computing the reflection function for any p, po, and cp. This function can then be 
used to calculate the reflected intensity for any directions of illumination and reflection and to 
find the plane, /fp(uo), and the spherical, A j, albedos: 


j i 2* i 

A P ( p 0 ) = — Jdpp Jd(pi?(p,p 0 ,cp) = 2 JV(p,p 0 )F dp. 


( 11 ) 


0 0 


A s = 2 j^(p 0 )p 0 dp 0 . (12) 

o 

Note that the BRF thus computed satisfies the fundamental principle of reciprocity: 

^(F.Fo’ ( P) = ^(Fo-F>9)- ( 13 ) 

Furthermore, since j(±p) = 1 if ra = 1, enforcing Eq. (8) ensures energy conservation for a semi- 
infinite nonabsorbing medium by rendering the plane and spherical albedos equal to 1 [cf. Eqs. 
(8), (II), and (12)]. 
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3. NUMERICAL ASPECTS AND COMPUTER CODES 


In this section we discuss numerical aspects and a FORTRAN implementation of the 
technique briefly outlined in the previous section. All computer procedures described are openly 
available on the World Wide Web at http://vvwvv.giss.nasa.gov/~-crmim/brt.html. 


3. I . Legendre expansion of the phase function 

The widely used Henyey-Greenstein phase function and its Legendre expansion coefticients 
are given by the following simple formulas: 


n&)-.. , * 2 , 3 , 2 . s<u ui, 

(l-2gcos0 + g ) 

(14) 

a s =(2s + \)g s . 

(15) 


Note that 


g =< cos© > 


(16) 


where 


1 + * l 

< cos© >= — \ / , (0)cos©d(cos0) = — a, 

2-i 3 


(17) 


is the asymmetry parameter of the phase function. Similarly, for the often used double-peaked 


Henyey-Greenstein phase function 1 

P(Q) = f 


60 


\ S\ + (1 - f) - — — ^ — 

(1 - 2 g, cos © + g, 2 ) 3/2 (1 - 2 g 2 cos© + g\ ) 3/2 


(18) 


with a positive g, and a negative g 2 , the Legendre expansion coefficients are given by 


=/ a *i+0-/X 2 ’ 


( 19 ) 


where a 5 , and a^ 2 are given by Eq. (15) with g = g, and g = g 2 , respectively. 
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The code for computing the expansion coefficients for polydisperse, homogeneous 
spherical particles is based on the standard Lorenz-Mie theory and the approach described in 


Refs. 61 and 62. The code allows one to select one of the following five size distributions: 


• the modified gamma distribution 


«(/•) = const x r a exp 


f ■/ \ 

ar' 

V Vc j 


( 20 ) 


• the log normal distribution 


n(r) = const xr exp 


f (Inr-lnr) 2 '] 


2 In' <7 


* ) 


( 21 ) 


the power law distribution 


n(r) = j const x r_3- r \ - r - r 2> 

[0, otherwise; 


( 22 ) 


the gamma distribution 


n{r) = const x r (1 36)/ *exp| 


f r_} 

cib 


, be (0,0.5); 


(23) 


• the modified power law distribution 



const, 0 < r < r,, 
const x(r/r,) a , r x 
0, r 2 < r. 




(24) 


The constant for each size distribution is chosen such that the size distribution satisfies the 
standard normalization 


J«(r)dr = l. (25) 

^iwn 
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Mathematically, particle radii in the modified gamma, log normal, and gamma 
distributions may extend to infinity. However, a finite /' max must be chosen in actual computer 

calculations. There are two practical interpretations of a truncated size distribution. First, r max 
can be increased until the scattering characteristics converge within some numerical accuracy 
(note that convergent r max values for the modified gamma and log normal distributions can be 

unrealistically large for small y or large c g 6j ). In this case the converged truncated size 

distribution is numerically equivalent to the distribution with r max = co . Second, a truncated 
distribution can be considered a specific distribution with scattering characteristics different from 
those for the distribution with r max =oo. Similar considerations apply to the parameter r min 
whose mathematical value for the modified gamma, log normal, and gamma distributions is zero, 
but in practice can be any number smaller than r max . Note that for the gamma distribution with 

r mjn = 0 and r max = oo, a and b coincide with the effective radius r elT and effective variance v eff , 

respectively, as defined by Hansen and Travis. 47 

The numerical integration of scattering characteristics over a size distribution is achieved 
by subdividing the entire interval [r mjn ,r max ] of particle radii into a number n of equal 

subintervals and applying a Gaussian quadrature formula with k division points to each 
subinterval. Note that n and/or k should be increased until the required numerical accuracy is 
reached. 

An efficient technique for computing the Legendre expansion coefficients for 
polydispersions of randomly oriented, homogeneous, rotationally symmetric nonspherical 
particles is described in detail in Ref. 64. This technique is based on the T-matrix approach 65 and 
an analytical method for averaging scattering characteristics over particle orientations. 66 
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The computation of the Legendre expansion coefficients for phase functions obtained 
with other numerical methods or measured experimentally is based on the numerical evaluation 
of the integral 

a s = f d0 sin QP(Q)P S (cos 0) , (26) 

- o 

which is a direct consequence of Eq. (6) and the orthogonality relation for Legendre 
polynomials. The integral is replaced by a Gaussian quadrature and an interpolation procedure is 
employed to find the phase function at Gaussian division points using the table of pre-computed 
or measured phase function values. The Legendre polynomials are computed using the 
recurrence relation and the initial conditions given by Eqs. (27) and (28) below with m = 0. We 
have found that spline interpolation usually provides quite acceptable results with the exception 
of phase functions having very sharp features, 67 such as the phase function for hexagonal ice 
crystals. The presence of the strong 22° and 46° halos in this latter case 68 necessitates the use of 
simple linear interpolation. Furthermore, the 5-function transmission peak in the ray tracing 
phase function for hexagonal ice crystals must be convolved with the Fraunhofer pattern, as 
described in Ref. 69. 

3. 2. Fourier components of the phase function 

The right-hand side of Eq. (5) is often written in terms of associated Legendre functions 

P s m (x) = (-i) m [(s + m)\/(s - m)\f /2 PQ m (x), i = V-T , rather than generalized spherical functions. 

It is well known, however, that the numerical computation of associated Legendre functions with 
large m and s is unstable and leads to overflows. 70 On the other hand, the computation of the 
generalized spherical functions via the upward recurrence relation^ 7,38 
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v ^7 + ir - m 2 C' (■'■) - (2s + ww - /?-#>“ /ffw 


(27) 


and initial conditions 


C;'(.r) = 0. CoW = (2i) 


/ (2m)! 

w!m! 


(1--V-/ 


(28) 


is numerically stable and efficient. Furthermore, the concept of generalized spherical functions 
naturally appears in the theory of polarized radiative transfer, 3 58 the Lorenz-Mie theory, 62 and 
the T-matrix method 66 and provides a natural and appealing link between these theories. 


3. 3. Iterative solution of the Ambartsumian 's equation 

By using a quadrature formula on the interval p e [0, 1] with n division points p^ and 

weights w , we convert integral equation (3) into a system of nxn nonlinear algebraic 
equations 


TU 


(P, +P*)« (P p »P,) = — ** (“Pp.P*) 




5=1 




(29) 


5 = 1 




j=i j'=i 


for the unknowns i? m (p ,p A p,q = . This system is solved by simple iterations using 




rn 




P m (-» p ,n q ) 


(30) 


as the initial approximation. The iterations converge very fast for m > 0 as well as for 
m- 0 and ra < 0.8 . However, the convergence rate becomes very slow when m = 0 and 1 - ra« 
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I. To accelerate convergence, we use a procedure similar to those developed in Reis. 51 and 52 
and based on the Sobolev -van de Hulst relation. Specifically, after each iteration, we compute 
the quantities [ef. Eq. (8)] 

Af , 1 (f l p ) = f (~l l p ) ~ > ( J ^) 

<,=i 

where j numbers iterations. We then improve /? [ ° 7 |(p / ,,(.i (/ ) by replacing it with 

/ C(*VM + K^jCPp)/^) + '(fip) A U](M]’ (32) 

where k is an appropriately chosen constant. This improved approximation is substituted in Eq. 
(31) to compute a new set of quantities A [yl (p p ), which are used again to further improve 


u„) via Eq. (32), and this procedure is repeated until 



(33) 


where e is the predefined absolute accuracy of computations. The improved y'th iteration 
^(y](Ep’. u ^) > s then substituted in the right-hand side of Eq. (29) to obtain the (/' +l)th 

approximation , which is again improved using Eqs. (31) and (32), and this entire 

process is repeated until 


max 

p.<7=i 




<€. 


(34) 


de Rooij 52 suggested to use the same value k = 0.5 in all cases. We have found, however, that 
this value may cause divergence when w < 0.995 and that k should be single-scattering-albedo 
specific. After having performed many numerical experiments, we have chosen the following k 
values: 
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K= < 


(35) 


0.5 for nr >0.995, 

0.1 for 0.95 < uj < 0.995, 

0.05 for 0.8 < ra < 0.95. 

The use of reciprocity [Eq. (13)] reduces the number of unknowns in Eq. (29) by a factor of 
2iV/(iV + 1) and provides a significant saving of computer resources. 

This numerical procedure renders only reflection Junction values /?(j.i. p .u (/ ,(p) at the 

division points of the quadrature formula. BRF values for p and po not coinciding with one ot 
the quadrature nodes must be found by numerical interpolation/extrapolation, which may result 
in lower accuracy than for the BRF values at the quadrature nodes. Therefore, the number of 
quadrature division points n should be increased until the desired numerical accuracy for all 
required BRF values is achieved. The accuracy can be significantly improved and n can be 
decreased by using the separation of the first-order scattering procedure (Section 3.7). 


3.4. Computation ofi(p) 


Taking into account the normalization 


2 


T l 

J P°(u,u')du' = 


(36) 


we derive from Eqs. (9) and (10) 

2(1 - ©) 


k =■ 


+i 

ra J u i{u) dz< 


(37) 


-l 


Replacing the integrals in Eqs. (9), (10), and (37) by respective quadrature sums, we obtain 


K±\i 0 ) = — - 


ns 


2(1 


X [‘(Vt/ ) p ° tePp- ' ^ ) + '(“Ik/ ) ■ p ° (iBp v )] . (38) 
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( 39 ) 


yZ v ‘>[ / ^) + / ^ l /^] =1 

“ />= i 

*=_ . ( 40 ) 

P =l 

Substituting /t = Vl-tn , /(p p ) = 2 , and /(-u p ) = 1/2 as the initial approximation, we calculate 
the right-hand side of Eq. (38) to obtain the next approximation for /(p p ) and /(-p p ) . Since 
this approximation may not satisfy the normalization of Eq. (39), we improve /(p p ) and /(-p p ) 
by dividing them by 

fi/'d'ag+'X-Ml- (4I) 

^ P = 1 

This improved approximation satisfies Eq. (39) and is used to compute the next approximation 
for k via Eq. (40). The new k, /(p p ) , and /(-p p ) values are substituted in the right-hand side of 

Eq. (38) to obtain the next approximation for /(u p ) and /(— u p ), and so on. The process is 
continued until /(p, ) and /(— p p ) converge within O.Ig. Note that this scheme is different from 
that described in Ref. 51. Dlugach and Yanovitskij 51 compute i(±p) using the expansion 

coefficients cx^ and a method of continued fractions. We have found, however, that the use of 
the expansion coefficients of the original phase function to compute /(±p) conflicts with the use 
of the renormalized phase function (Section 3.6 below) in Eq. (29) and may lead to divergence of 
the iterative solution of Eq. (29) for highly anisotropic phase functions. Our new procedure for 
computing /(±p) uses the already renormalized phase function and produces numerically stable 

and convergent results. 
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J. 5. Numerical integration 

The Gauss quadrature formula (e.g.. Ref. 71) has the highest algebraie degree ot precision 
(i.e., a formula with n nodes is exact for all polynomials of degree In - 1 and lower) and is 
traditionally used in radiative transfer calculations to numerically evaluate integrals on the 
interval [0.1] (e.g.. Refs. 72-74). A significant disadvantage of this quadrature is that the largest 
node is always smaller than 1, and if BRF values tor normal incidence and/or reflection are 
required, then one must use an extrapolation procedure. Unfortunately, extrapolation often 
produces poor numerical accuracy (e.g., see discussion on pages 210 and 211 of Ref. 52) and 
necessitates the use of the Gaussian quadrature formula with a large number of nodes. We have 
found that a more efficient approach is to use the so-called Markov quadrature formula (Chapter 
9.2 of Ref. 71) with one predefined node at p = 1. This formula still has the highest possible 
algebraic precision and is exact for all polynomials of degree 2n - 2 and low r er. Furthermore, it 
allows one to avoid the use of the extrapolation procedure or so-called extra points (Ref. 47 and 
L. D. Travis, personal communication). Multiple numerical tests have shown that the Gaussian 
and Markov quadratures with a number of nodes n larger than about 10 have essentially the same 
numerical accuracy for intermediate p and po values, whereas the Markov quadrature produces 
much better accuracy for p and/or p 0 equal to 1. Since the CPU time consumption in solving Eq. 
(29) is proportional to the use of the Markov quadrature with a reduced number of nodes 
results in a significant saving of computer resources. We have developed a simple, efficient, 
and highly accurate FORTRAN procedure for computing the nodes and weights of the Markov 
quadrature formula with an arbitrary n. Table 1 exemplifies the performance of the procedure 
and lists the nodes and weights of the Markov quadrature with n = 30. 
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The direct application of a quadrature formula to the integration ^-interval [0. 1] is a 
standard approach in the radiative transfer theory (e.g., Rets. 47-49, 51, 55, 61, 72-74). 
However, it provides poor sampling of zenith angles close to 0° and, as multiple numerical tests 

have shown, causes a very slow convergence of /?°(l,l) with increasing n for particles large 
compared to the wavelength. This happens even when the Markov quadrature is used and no 
extrapolation is involved or when the Gaussian quadrature is used along with an extra point at 
p = ji 0 = 1 . On the other hand, convergence with increasing n is fast for p < 1 and p 0 < 1 . We 

have found that a very efficient way of avoiding excessive n values in radiative transfer 
computations is to apply the Gaussian quadrature to the interval [0, 7i/2] of zenith angle values. 
Since 


71/2 


| /(p)dp = | /(cos§)sin$d9, 


(42) 


we easily derive the following expressions for the respective division points and weights in Eq. 
(29): 


[ 71 71 

^= cos 4^+4 


\v p = W p sin 


7t __ 7C « 

-X +- , p = \,...,n , 
4 p 4 J 


(43) 


where X p and W p are Gaussian nodes and weights, respectively, on the interval [-1,+ 1] . These 
division points provide a much better sampling of zenith angles close to 0° and a much higher 
convergence rate for i?°(l,l) with increasing n than the Gaussian or the Markov quadrature 
formula applied to the interval [0, 1] of p values. 
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3.6. Renormalization of the phase function 

Although analytically the zeroth Fourier component of the phase function must be 
normalized according to Eq. (36), the numerical evaluation of the left-hand side of Eq. (36) 
usually produces p-dependent numbers not equal to 1: 


2 


X 

t?=l 


\V 


</L 


/>V 


.pj + p (p„ -p. 



* 1 . 


(44) 


This results in a deviation of the “numerical - ' single scattering albedo from its actual value and, 
for nonabsorbing or weakly absorbing media, can lead to an efficient “photon gain” or “photon 
loss.” A direct adverse consequence is a serious violation of energy conservation and poor 
numerical accuracy. Hansen 55 developed a so-called renormalization procedure, which 

numerically enforces the normalization of Eq. (36) by slightly modifying the / ,0 (p / ,,p ? ) values. 


We have found that the renormalization procedure of Ref. 55 produces accurate BRFs in most 
cases, but not always. Therefore, we have developed an alternative renormalization procedure, 
which is simpler than that of Ref. 55 and appears to be more stable. Specifically, we multiply 

the quantities P°(a p ,\i p ). p = by the correction factors 


£ =1 + 


2-25, 


w p P (p^pJ 




(45) 


This correction makes the left-hand side of Eq. (44) equal to 1 for any p and is applied to higher 
Fourier components of the phase function as well. Since it affects only the forward-scattering 
values of the phase function, it has negligible effect on the bidirectional reflection function while 
numerically ensuring energy conservation. 
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3. Separation of the first-order-scattering contribution to the refection function 

For large scattering particles with highly variable phase functions and for p and p 0 
significantly smaller than 1, one may need very many Fourier terms in Eq. (2) in order to 
accurately represent the reflection function. On the other hand, it is also known that with p and 
j-io approaching zero the only significant contribution to the reflection function comes from 
photons scattered only once .' 6 This suggests the idea of subtracting the first-order-scattering 
contribution from all Fourier components of the reflection function, thereby greatly reducing the 
number of numerically significant Fourier components, evaluating the right-hand side of Eq. (2), 
interpolating (if necessary) this slowly varying high-order-scattering part of the reflection 
function, and finally adding the exact single-scattering contribution . 70 The latter contribution can 
be easily computed analytically for the scattering angle 0 corresponding to a specific 
combination of p, po, and cp values and given by 

cos© = -pp 0 + Jl -p 2 -Jl -p 2 cosep. (46) 

In other words, the total reflection function is represented in the form 

/?(p, p 0 , 9 ) = R { (p, P 0 , <p) + X (2 “ 6 m 0 )[tf W (P>Po) - K (IT P 0 )] cos m S> , (47) 

m - 0 

where 5 mm ’ is the Kronecker delta, 


^i(ITPo’ ( P)= P(®), 

4(P + P 0 ) 

(48) 

*r(n,m>)- ,, o)> 

4(p + p 0 ) 

(49) 
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and /;/, «/w nm . 48 The term in square brackets on the right-hand side of Eq. (47) is a smooth 
function of p and po and can be accurately interpolated even when the number ot quadrature 
nodes is relatively small, while P( 0) is computed via Eqs. (6) and (46). 


4. COMPUTATIONS AND DISCUSSION 

4. 1. Soil surfaces 

Table 2 lists parameters of four soil particle models used in the computations described 
below. We assume the standard gamma size distribution of Eq. (23) with an effective radius of 
a = /- e|T = 10pm and an effective variance of b = v eff = 0.1 . This effective radius is typical of soil 

particles (e.g.. Ref. 76). The four values of the refractive index m = m r +irn t are also typical of 

soil particles at the visible wavelength X - 0.63 pm considered. 76 The single scattering properties 
of the four soil particle models were computed assuming the spherical particle shape and using 
the Lorenz-Mie theory. (It should be noted, however, that the Lorenz-Mie theory does not 
necessarily provide the best representation of soil particle phase functions. 77 so ) Table 2 gives the 
respective values of the single-scattering albedo vs , asymmetry parameter of the phase function 
< cos© > , the number of terms in the Legendre decomposition of the phase function s max [Eq. 

(6)], and the spherical albedo A s . Note the significant decrease of gt and increase of < cos© > 

with increasing imaginary part of the refractive index. The solid curves in the upper panel of 
Fig. 2 show the respective Lorenz-Mie phase functions, while the dotted curves show the 
asymmetry-parameter-equivalent Henyey-Greenstein phase functions [Eq. (14)]. 

Table 2 also gives spherical albedo values computed using the equivalent Henyey- 
Greenstein phase function, /^(HG) , and the simple approximate formula 
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.-Jy(HH) = (1 - v)/(l + .v) 
derived by Hovenier and Mage, 81 where 


( 50 ) 


5 = 

is the so-called similarity parameter.' 6 It is seen that the My(HG) values are quite close to the 
exact ones, while the Hovenier and Hage approximation provides somewhat lower accuracy. 

Solid curves in the upper panel of Fig. 3 depict the plane albedo A P as a function of the 
cosine of the illumination zenith angle ji 0 . Note that A p is determined only by the 0th 

component of the reflection function via Eq. (11) and, as a consequence, the computation of the 
upper panel of Fig. 3 using 50 quadrature division points took less than 2 sec of CPU time on an 
IBM RISC model 397 workstation. We also computed the plane albedo using the asymmetry- 
parameter-equivalent Henyey-Greenstein phase functions and the 5-Eddington approximation. 
The ratios of these approximate plane albedo values relative to the exact ones are shown by the 
dotted and solid curves, respectively, in the bottom panel of Fig. 3. Not surprisingly, plane 
albedos decrease with increasing the imaginary part of the refractive index and, thus, decreasing 
the single-scattering albedo. Both the 8-Eddington approximation and the asymmetry-parameter- 
equivalent Henyey-Greenstein phase function produce significant errors, especially for grazing 
illumination. Using the asymmetry-parameter-equivalent Henyey-Greenstein phase function 
overestimates the plane albedo for small p 0 and underestimates it for p 0 close to 1, which is 

naturally explained by the scattering-angle pattern of the phase function differences (upper panel 
of Fig. 2). The errors increase significantly with increasing absorption. This can be explained 
by the increasing contribution of photons scattered only once and by the large differences in the 
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single-scattering phase functions. The much better accuracy of the .-Jy(IlG) values in Table 2 

can be explained by cancellation of the plane albedo errors after integrating over p 0 in Eq. (12). 

Figure 4 shows the angular distribution of the reflected intensity computed lor soil 
particle model 1 using the exact Lorenz-Mie phase function and its Henyey-Greenstein 
counterpart and assuming F= 1 in Eq. (1). The computations tor the exact phase function used 
100 quadrature nodes and took about 25 min of CPU time including the Lorenz-Mie computation 
of the Legendre expansion coefficients, the solution of Eq. (3), and interpolation. The 
computations for the equivalent Henyey-Greenstein phase function used 50 quadrature nodes and 
took less than 3 min. 

Two obvious features of the reflected intensity distributions shown in the left column are 
the backscattering enhancement (p = p 0 ,cp = 180°) caused by the glory in the Lorenz-Mie phase 
function (upper panel of Fig. 2) and the strong near-forward scattering for the cases of grazing 
and near-grazing incidence caused by the phase function peak at small scattering angles. The 
reflectance patterns for the equivalent Henyey-Greenstein phase function exhibit only the second 
feature, which is explained by the absence of the backscattering phase function peak similar to 
the glory. The right column of Fig. 4 shows that BRF errors caused by the use of the equivalent 
Henyey-Greenstein phase function can be very large and can exceed a factor of 20 at 
backscattering geometries and a factor of 3 at near-forward-scattering geometries. These errors 
can be unequivocally attributed to the phase function differences. Thus, Fig. 4 makes a strong 
case against using approximate phase functions in BRF computations for semi-infinite 
particulate media. 
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4.2. Snow surfaces 

In this section we describe BRF computations for three snow particle models summarized in 
fable 3. The assumed wavelength is X = 0.65 pm. Model 1 particles have highly irregular, 
random-fractal shapes introduced in Ref. 83; model 2 particles are homogeneous ice spheres; and 
model 3 particles are regular hexagonal ice crystals with a length-to-diameter ratio of 2. The 
nonspherical model l and 3 particles are randomly oriented in three-dimensional space. For all 
three models we used the same power law distribution of radii or projected-area-equivalent- 
sphere radii [Eq. (22)] with an effective radius of 50 pm and an effective variance of 0.2. The 
respective phase functions were computed using the ray tracing technique 8j> coupled with the 
Kirchhoff approximation 69 for models 1 and 3 and the Lorenz-Mie theory for model 2. They are 
shown in the lower panel of Fig. 2 and exhibit large differences exceeding an order of magnitude 
at some scattering angles. The resulting differences in the asymmetry parameter are also 
significant (Table 3). As discussed in Ref. 84, the phase functions of models 1 and 3 may 
represent limiting cases of highly distorted and pristine ice crystals, respectively. Since water ice 
is essentially nonabsorbing at visible wavelengths, the single-scattering, plane, and spherical 
albedos for all three models are equal to unity. 

Figure 5 shows the reflected intensities for the three snow particle models, while Figure 6 
depicts the ratios 2/1, 3/1, and 3/2 of intensities for the respective models. It is obvious that the 
shape of the scattering particles has a profound effect on the reflectance of flat snow surfaces. 
Although radiance differences between the different models are relatively small at nearly normal 
incidence (po = 0.9), they significantly increase with decreasing po and result in intensity ratios 
smaller than 0.2 and larger than 3 (cf. Ref. 84). This is a direct consequence of the increasing 
relative contribution of the first-order scattering to the reflection function and the large phase 
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function differences. Hexagonal snow crystals (model 3) produce the most structured radiance 
field dominated by the baekscattering peak and the primary (22°) and secondary (46°) halos (two 
lower panels in the right column of Fig. 5). These features clearly show up in the 3/1 and 3/2 
intensity ratios (Fig. 6). The spherical ice particles produce a noticeable enhancement ot 
intensity caused by the rainbow. This feature is especially well seen in the 2/1 ratio. The 
radiance field produced by the featureless phase function of irregular snow crystals (model 1) is 
by far the least structured (left column of Fig. 5). These results emphasize the importance of 
accurate treatment of single-scattering phase functions for realistic snow grain models. 

5. CONCLUSIONS 

We have described in detail an efficient technique for computing bidirectional reflectance 
of semi-infinite discrete random media based on an exact numerical solution of the radiative 
transfer equation. This technique results in a very compact and fast computer code and produces 
BRFs which fully comply with reciprocity and energy conservation. The high efficiency and 
accuracy of the technique make far less tempting the use of approximations such as the 5- 
Eddington approximation, the asymmetry-parameter-equivalent Henyey-Greenstein phase 
function, and the truncation of the phase function 73 and provide an ideal tool for testing the 
accuracy of approximate approaches. 81,85 Our sample computations for flat soil and snow 
surfaces have clearly demonstrated the limited applicability of approximate treatments of the 
single-scattering phase function in BRF modeling. 

Since we considered only nonabsorbing or weakly absorbing media, we ignored the 
opposition effect caused by so-called shadow hiding. 1 Other factors ignored in the model are the 
effects of polarization, 53 ' 33 packing density, 33 ' 41 ’ 86-91 and coherent baekscattering. 29,30 ' 32,38 ’ 43,92 
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However, they can be incorporated in a rather straightforward manner, as described, e.g., in Refs. 
26-28, 52, and 93-95. This is the subject of our current research. An equally challenging 
problem of macroscopic surface roughness 812 ’ 19,24 ’ 46 ’ 96 can be addressed by convolving BRFs 
computed as described in this paper with a Monte Carlo procedure handling multiple surface 
reflections and surface shadowing. As pointed out in the introduction, this approach avoids 
time-consuming ray tracing inside a nonabsorbing or weakly absorbing, optically semi-infinite 
medium and provides a great saving of computer resources. The ultimate challenge is to take 
into account the effects of the discontinuous nature of light scattering in densely packed discrete 
random media, 31 ' 35 97 but this requires the development of a much more sophisticated approach. 
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FIGURE CAPTION 


Fig. I. Spherical coordinates specifying the direction of light propagation. The zenith and 
azimuth angles of the incident beam are O 0 > tc/2 and <p 0 = 0, respectively. 

Fig. 2. Phase functions for soil particle models 1—4 (upper panel) and snow particle models 1-3. 
Dotted curves in the upper panel show asymmetry-parameter-equivalent Henyey-Greenstein 
phase functions. 

Fig. 3. Upper panel: plane albedo versus go for soil particle models 1-4. Lower panel: plane 
albedos computed using the 5-Eddington approximation (solid curves) and asymmetry- 
parameter-equivalent Henyey-Greenstein phase functions (dotted curves) relative to the exact 
values for soil particle models 1-4. 

Fig. 4. Left column: reflected intensity versus p and cp for soil particle model 1. Middle 
column: the same but for the asymmetry-parameter-equivalent Henyey-Greenstein phase 

function. Right column: the ratio of the intensities shown in the middle and left columns. The 
four values of the cosine of the illumination zenith angle po = 0.9, 0.7, 0.4, and 0.1 are indicated 
by the yellow stars in the right column. The azimuth angle of the incident radiation is zero. 

Fig. 5. Reflected intensity versus p and (p for snow particle models 1, 2, and 3. The four values 
of the cosine of the illumination zenith angle po = 0.9, 0.7, 0.4, and 0.1 are indicated by the 
yellow stars. The azimuth angle of the incident radiation is zero. 

Fig. 6. As in Plate 4, but for ratios of reflected intensities. 
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Table 1 . Division points and weights of the Markov quadrature 
formula on the interval [0. I ] with n - 30 




1 0.00160587785254 0.00411899413797 

2 0.00844194880403 0.00954444032329 

3 0.02066193363616 0.01487343889309 

4 0.03813469699711 0.02004018939073 

5 0.06066914207658 0.02498751843953 

6 0.08801845813989 0.02966112217483 

7 0.11988302768490 0.03400976899322 

8 0.15591374708438 0.03798580509944 

9 0.19571586157012 0.04154566451353 

10 0.23885329363337 0.04465034297146 

11 0.28485342210352 0.04726582410455 

12 0.33321226086401 0.04936345167913 

13 0.38339998090611 0.05092024336491 

14 0.43486671538603 0.05191914244132 

15 0.48704858415304 0.05234920462588 

16 0.53937387177272 0.05220571795295 

17 0.59126929137078 0.05149025437772 

18 0.64216626567319 0.05021065253237 

19 0.69150715642309 0.04838093181641 

20 0.73875137391194 0.04602113875159 

21 0.78338129966059 0.04315712726339 

22 0.82490795730055 0.03982027524838 

23 0.86287636938624 0.03604714041 164 

24 0.89687054110030 0.03187905880664 

25 0.92651801527194 0.02736168939988 

26 0.95149394561648 0.02254450560005 

27 0.97152463085599 0.01748022010489 

28 0.98639040253929 0.01222402012888 

29 0.99592721636071 0.00683100534121 

30 1 .00000000000000 0.00 1111 11111111 
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Table 2. Soil particle models 


Model 

a (mp) 

b 

m r 

m, 

w 

< COS0 > 

^ max 

-‘s 

A,{UG) 

/Is(HH) 

1 

10 

0.1 

1.55 

0.001 

0.85404 

0.83752 

641 

0.1399 

0.1382 

0.1655 

i 

10 

0.1 

1.55 

0.002 

0.76137 

0.86568 

644 

0.0727 

0.0716 

0.0889 

3 

10 

0.1 

1.55 

0.003 

0.69923 

0.88582 

645 

0.0472 

0.0464 

0.0588 

4 

10 

0.1 

1.55 

0.004 

0.65646 

0.90054 

646 

0.0345 

0.0339 

0.0435 


Table 3. Snow particle models 


Model 

Shape r 

iff C um ) 

V eff 

m r 

m, 

w 

< COS0 > 

^max 

1 

Irregular 

50 

0.2 

1.311 

0 

1 

0.7524 

2000 

2 

Spherical 

50 

0.2 

1.311 

0 

1 

0.8860 

1948 

3 

Hexagonal 

50 

0.2 

1.311 

0 

1 

0.8117 

2000 
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Plane Albedo 























